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ABSTRACT 

The primordial elements, H, He and Li are included in a low temperature equation of 
state and monochromatic opacity function. The equation of state and opacity function 
are incorporated into the stellar evolution code NG-elms, which makes use of a non- 
grey model atmosphere computed at runtime. NG-elms is used to compute stellar 
evolution models for primordial and lithium free element mixtures, for stars in the 
sub-Solar mass range 0.8-0.15 Mq. We find that lithium has little or no effect upon 
the structure and observable properties of of stars in this mass range. Furthermore 
lithium is completely destroyed by fusion before the main sequence in stars of mass 
less than ~0.7 Mq. We find that on the red giant branch and Hayashi track, the use 
of a non-grey model atmosphere to provide the upper boundary conditions for the 
stellar evolution calculation, results in significantly cooler less luminous stars, across 
the mass range. 
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1 INTRODUCTION. 

, The observation and identification of the first generation 
■ of stars (Population III) re mains one of the main go als of 
I modern astronomy. Indeed, iKashlinskv et all i2004ft have 
investigated the possibility that light from the first stars 
may contribute to the near infrared background. They also 
suggest that this light may h ave already been observed 
via the Spitzer space telescope iKashlinsk v 2006). However 
iMadau fc Sil^ ll2005r cast doubt upon this conclusion. 

If population III stars formed with masses less than 
0.8 M© they would still be on the main sequence at the 
present day. This leads to the possibility that these Popula- 
tion III stars may be present and observa ble within the so- 
lar ne ighbourhood. IChristlieb etal] ll2002h and lFrebel et alJ 
i2005ll have reported the discovery of two low mass stars 
with [Fe/H] < -5. It has been suggested that these stars 
may be population III stars which have accreated metals ei- 
ther from the interstellar medium of from a companion star 
to achieve their present day metallicities. 

In the absence of metals most current star for- 
mation models predict that the Population III ini- 
tial mass functio n is skewed to the ve r y hig h mass 
(> 100 MfT)) jBromm. Coppi fc Larsonl 119991 l2002t 
Abel et alJl2002h. However this scenario is not yet certain, 
Tsuiimoto fc Shigevamal (|20o3) suggest that intermediate 
mass (3.5 - 5 Mq) primordial stars are required to produce 
the heavy element abundance ratios seen in some metal 



poor stars. Furthermore, iNaeakura fc Qmukai (|200^ have 
proposed a method by which low mass primordial stars 
can form from the unpolluted remnants of high mass star 
formation, if the high mass star forms a black hole without 

supernova. 

In our earlier paper jHarris et alJ [iooii), we com- 
puted evolutionary models of low-mass metal-free stars, us- 
ing only a hydrogen and helium equation of state. It was 
demonstrated that HJ has a strong effect upon the ef- 
fective temperature and luminosity of stars with M<0.4 
M0. However, these models were computed using grey 
model atmospheres to provide the outer boundary condi- 
tion for the equations of stellar structure. For very low 
masses (M <0.5 Mq) the grey approximation is known to be 
inaccurate feaumo n et al ]|l994l: iBaraffe fc Chabrieil Il996l : 
pJaraffc ct al. 1997: Chabrier fc Baraffelll997^. Additio nally, 
our models neglected lithium, ' Maver fc DuschJ |2005") have 
shown that Li can act as a significant electron donor at low 
temperatures in a primordial (H, He and Li) mixture. In this 
work we update our equation of state to include Li, and com- 
pute evolution models with a non-grey model atmosphere, 
from protostellar collapse to the Helium flash. 



2 NON-GREY EVOLUTION. 

The non-grey evolution of low mass stars (ng-elms) code is 
based upon a combination of 2 existing codes, CESAM iMoreJ 
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Il997l) and marcs iGustafsson et aljll975l) . cesam is a one 
dimensional stellar evolution code developed for the mod- 
elling of low and intermediate mass stars, which makes use 
of a grey model atmosphere. MARCS is a plane parallel model 
atmosphere code, for which we have rewritten the opacity 
and equation of state subroutines and modified to run as a 
subroutine of cesam. marcs provides the surface boundary 
conditions for CESAM at runtime. This code will be reported 
in more detail in a future article. 

To provide nuc lear reaction rates we use the NACRE 
jAneulo et alJll999ll compilation, for the PP chain, CNO 
cycle and triple alpha process. The following reactions are 
accounted for: H(p,e+ i/)D, D(p,7)3He, 3He(3He,2p)*He, 
^He(3He,7)''Be, ■^Li(p,a)*He, ■'Be(e-,i^ 7)'^ Li, 
^Be(p,7)»B(e+ z/)*Be(a)^He, ^'^C{vaf''^{e+ vf'C, 
i3c(p,7)"N, "N(p,7)i5o(e+ vf^'^, i5j^(p,7)i6o, 
i^N(p,a)i2C, ^«0(a,7)"F(e+ vfO, ^'0(p,a)"N, 
''m{2aaY^C, i2C(q,7)1''0, ^''0(Q,7)2''Ne. 

At low temperatures (T<9000 K) we use the equation of 
state and opacity function discussed in sections |2 . II and 12.21 
A linear interpolation upon the the OPA L equation of state 
jRo gers fc Navfonovl |2003) and opacity (Iglcsias & Rogers 
Il996f) tables provides the equation of state and opacity data 
at higher temperatures (T > 11 000 K). Between tempera- 
tures of 9000 and 11 000 K a linear interpolation between 
OPAL data and the data of this work is used. The OPAL 
equation of state tables, cover the temperatures and densi- 
ties of stars with masses greater than about 0.15 M0. This 
limits our calculation of evolution models to masses of 0.15 
M0 and greater. Conductive opacities arc obtained by inter- 
polating the tables oflHubbard & Lampc ( 1969). The diffu- 
sion of ele ments including met a ls is a ccounted for via the 
method of iMichaud fc ProfBtd (|l99.'t . The mixing length 
theory of convection is used, with an adopted value of the 
mixing length parameter of a = 1.745. This value of the 
mixing length parameter was achieved via solar calibration, 
using an equation of state which contains the elements H, 
He, Li, C, N, O, Ne, Na, Mg, Al, Si, S, Ca, Ti, Fe; this wiU be 
fully reported in a future article. At the accepted solar age 
of 4.5 Gyr, for a solar mass model we achieve Teff = 5781.8, 
L/Lq = 1.00398 and R/Rq = 1.000497. Our initial mass 
fractions are Xq = 0.70440, Zo = 0.018698, and surface 
mass fractions at the solar age are Xs = 0.72889, Zs = 
0.017783, giving Z^jX^ = 0.02440, t his is comparable wit h 
the solar value of Zq/Xq = 0.02440 l|Grevesse et alJll99fi() . 
These initial conditions compare with the solar evolution 
and h elioseismology model fits of iGuenther fc Demarauel 
who find Zo=0.020 and Xq = 0.7059, and achieve 
Zq/Xq =0.0244. 

During the cour se of this work 

lAsplund. Grevesse fc SaxivaJ l^2005^ have published new 
estimates of the metal abundance ratios and metallicity of 
the Sun. They revise downward the metal to hydrogen mass 
fraction ratio to Z/X = 0.0165. As metals are important 
contributors to the solar opacity, the reduced metallicity 
will result in a lower opacity. So that a solar model com- 
puted with the new metal mix will have a greater effective 
temperature, luminosity and radius than o ne computed 
with the metal mix of jGrevesse et a l .1 11 9961) . To counter 
this, the mixing length parameter will need to be revised 
downward. 



2.1 Low temperature equation of state. 

In order to compute opacities it is first necessary to know 
the abundances of each atomic molecular and ionic species. 
The equation of state computes the various thermodynamic 
quantities as well as the number densities of each species. 
Our low temperature equation of state (lomes) is an ex- 
tensively rewritten version of the equation of state used in 
iHarris et alJ ll2004afl. This equation of states has been used 
by 'Eneel et al' ('2005") and in a modified, non-ideal, form 
by Harris et al. (2004b) to demonstrate the importance of 
HeH""" at very low hydrogen ab undances. In light of the find- 
ings of lMaver fc DuschJ l^2005^ . we have updated the equa- 
tion of state to include Li, Li^, Li""""^, LiH and LiH"^. lomes 
is based upon the Saha equation. Each species is written 
in terms of its constituent neutral atomic species and elec- 
trons. For example the formation reaction of H^ is written 
as SH^Hg'-l-e", the Saha equation is then 

N{Ht)N{e-) _ QT{H+)QT{e-) ( AE\ 

[N{H)]--^ [QriH)]-^ ^\ kTj ^ ^ 

where A'' are the number densities of each species, Qt are 
the total partition functions, AE is the energy difference 
between the left and right hand sides of the reaction equa- 
tion. The Saha equation for each species is used to con- 
struct simultaneous equations for the conservation of charge 
and conservation of nucleons, which are th en solved. This 
is a widely used tech n ique, see for example iKurucd lll97Cfl 
and iMaver fc Duschll ll2005ll . The rewritten a lgorithm has 
none o f the convergence problems reported in lHarris et alJ 
(Hooii). The primordial atomic, molecular and ionic species 
included in the equation of state are quoted in tabled For 
the mass range covered here, the densities of the regions 
below 9000 K are not sufficiently high to result in strong 
non-ideal efi^ects such as pressure dissociation and pressure 
ionisation. Thus a Saha based equation of state is adequate 
for the purposes of this study. 

Figure shows how number density of the dominant 
positive ions and electrons change with temperature for a 
density of 10"'' g cm~^ . As found bv lMaver fc DuschJ ll2005l) 
ionised lithium is the dominant positive ion at low tempera- 
tures, above 3000 K the abundance of H J exceeds that of Li^ 
and H""" becomes the dominant positive ion for temperatures 
greater than 3600 K. At densities lower than lO"'^ g cm"'', 
the abundance of H^ is hugely reduced by free electrons 
from Li, to the extent at which it becomes unimportant. At 
these low densities, Li""" is the dominant positive ion at low 
temperatures and H^ is dominant at high temperatures. 

2.2 Low temperature opacities 

Our sources of continuous opaci ty and s c atterin g data are 
the same as those used in Ha rris et al ] ^M4^. Most of 
our computations involve non-grey atmospheres which re- 
quire monochromatic opacities over a range of frequencies. 
The use of line opacity i n the comput ation of model at- 
mospheres is problematic ^ICarbonlll974^ ■ by its nature line 
opacity is not smooth. Spectroscopic lines block electromag- 
netic fiux from escaping an atmosphere near the line centres, 
but allow fiux to escape between lines. Much work has gone 
into statistical treatments of line opacity within model at- 
mospheres, which has resulted in some success for opacity 
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Table 1. The primordial atomic, molecular and ionic species included in the equation of state and the source 
of the partition function used. 
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Figure 1. The number densities of the dominant positive ions and electrons at a density of 10~® g cm~ 



sampling an d the opacity distrib ution function (ODF), see 
for example lEkbere et al.l l)l98f^) . In ng-elms it is neces- 
sary for model atmospheres to be computed approximately 
30 times per evolutionary time step, it is therefore essential 
for a solution to be arrived at rapidly. With this in mind 
we use a straight mean opacity computed over the inter- 
val betvyeen frequ ency points, for each source of line opac- 
ity. ICarbonI il974ft has shown that the straight mean tends 
to over estimate opacity, we find that the systematic er- 
rors introduced by a straight mean are relatively small from 
many closely packed weak lines, such as molecular bands. 
However, for strong isolated lines, such as atomic lines, the 
opacity is dramatically over predicted. For the mass range 
studied here the atomic (H, Li) and molecular (H2) lines 
make up a relatively small fraction of the total opacity, we 
have therefore neglected H2, Li and H lines from our cal- 
culations. To test if line opacity from H2, Li and H, affects 
the evolution of low mass stars we have computed O DFs for 
these species. Line lists for H and Li are taken from iKurucd 
il99 3). Einstein A coefficients for H2 quadrupole transitions 
taken from Wolniowicz ct aT] lll998l) and H2 energy levels 
are taken from the Molecular Opacity Database UGAMOP 
l |http: / / www.p hysast .uga.edu / ugamop / I. Two evolutionary 



models are computed using these ODFs and are discussed 
in section |3] 

Figure 121 shows the continuous opacity and Rosseland 
mean opacity at 5000 K for densities of 2.5xl0~* and 
5.0 X 10""^ g cm~^. For the higher density (5x10"'' g cm~^) 
opacity is completely dominated by H~ bound-free and free- 
free continuous opacity. The continuous opacity is slightly 
less than the Rosseland mean for frequencies above 22 000 
cm~^, and slightly greater than the Rosseland mean for 
frequencies between 8000 and 22 000 cm~^. Thus, overall 
the Rosseland mean opacity, which is a weighted harmonic 
mean, gives a reasonable representation of opacity over the 
frequency range of the peak of the Planck function. At the 
lower density (2.5xl0~* g cm~'^), the opacity is still domi- 
nated by H~ in the red, but in the blue the contribution to 
opacity from Rayleigh scattering and H bound-free contin- 
uous opacity are important. The Balmer jump is noticeable 
and the Rosseland mean opacity appears to be too low to 
provide a good representation of overall opacity. As is dis- 
cussed in section 1231 this 'overestimate' of opacity leads to 
the breakdown of the grey approximation for giant stars. 

Radiative transport in the interior of the star is treated 
using the diffusion approximation, which requires the Rosse- 
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Figure 2. Continuous opacity at 5000 K, for densities of 2.5xl0~* and 5.0 X 10~® g cm"-^, which are the approximate densities found 
in giant and dwarf photosphere respectively. For reference, the Rosseland mean opacity is shown as a horizontal line for each of the two 
curves. 



land mean opacity. Our grey model atmospheres also use 
the Rosseland mean opacity. The Rosseland mean opacity 
is computed at runtime, via numerical integration using 145 
frequency points. 

2.3 Grey and non-grey model atmospheres. 

It has previously been found that for dwarf stars with 
Tcff less than ~5000 K, the grey approximation pro- 
vides inaccurate boundary condition s for stellar evolu- 
tion calculations (Saumon ct al."l994"; 'Baraffe fc ChabrieJ 
[l996; Baraffc ct al. 1997; Chabricr & Baraffc 1997,). Figure 
13 shows the temperature-pressure curves of 4 grey and non- 
grey model atmospheres, computed for Tcff=5000 K and 
with surface gravities of logg — 1.5 (giant) and logg = 5 
(dwarf). It is clear that for both surface gravities the grey 
approximation provides a poor representation of the non- 
grey atmosphere. In general the non-grey atmosphere is less 
dense than the grey atmosphere, this is consistent with the 
grey atmosphere having an effectively lower opacity than 
the non-grey atmosphere. There is a far stronger deviation 
of grey atmosphere from non-grey atmosphere at low surface 
gravities. 

As temperature is reduced the non-grey effect strength- 
ens. For TcfT above ~6000 K the grey atmosphere gives good 
agreement with the non-grey atmosphere, even for low sur- 
face gravities. Non-grey model atmospheres must be used to 
provide the boundary conditions for stellar evolution models 
below 6000 K, particularly for giants. 

2.4 Computation of evolution models. 

In the computation of evolution models we use two separate 
element mixtures. A lithium free mixture which contains 
only hydrogen and helium and a primordial mixture which 
in addition to H and He also contains Li. We adopt the 
element abundance ratios determined by ICoc et al] j20o3) 
from big bang nucleosynthysis model fits to the cosmic 
microwave background observations by the Wilkinson mi- 
crowave anisotropy probe (WMAP). It must be noted that 
th e value of the Li/ H ratio of 4.25lg4g ^ 10"^", determined 
bv lCoc et alJ ll2004l) , deviates by a factor of 3 from the value 



of 1.2 3_n X 10~^ d etermined from observation of halo 
stars (iRvan et al.ll200ol) . We adopt the higher value in order 
to asses the maximum impact of Li on low mass stellar evolu- 
tion in the early Universe. The adopted element abundances 
for the primordial mix is X = 0.75210 and Zli = 2.24x 10"^, 
and for the lithium free mix is X = 0.75210 and Zli = 0. 

All Models are computed from an initial quasi-static 
protostellar model undergoing gravitational collapse. As in 
our earlier work |h arris ot al. 2004a) we use a contrac- 
tion constant of 0.015 LqMq^K"^. The models are followed 
through collapse down the Hayashi track, through the main 
sequence. Models are allowed to evolve onto and up the red 
giant branch to near the point of helium flash, or the point 
at which stars of insufficient mass to fuse Li turn off the red 
giant branch. With ng-elms we have computed evolution- 
ary models for stars with masses of between 0.15 and 0.8 
M0, using grey and non-grey model atmospheres, for both 
primordial and Li free mixtures. 



3 RESULTS. 

Evolutionary tracks on a HR diagram, for stars of, 0.15, 
0.4, and 0.6 M©, computed using grey and non-grey model 
atmospheres, are shown in fig ure |1] As we have previ- 
ously found iHarris et al.ll200al . on the red giant branch 
and Hayashi track, all the non-grey models have signif- 
icantly lower effective temperatures for a given luminos- 
ity than their counterparts computed with grey model at- 
mospheres. Furthermore, on the Main Sequence the non- 
grey models with masses less than ~0.5 M0, are also 
cooler than thei r grey counterparts, w hich is consistent with 
the findings of 'Saumon et al.' ('1994'), hereafter SBLHB94, 
and of Bar affc fc C habri er (1996); Baraffe e t al. (1997|); 
IChabrier fc Baraffj lll997l) . These two factors again, illus- 
trate the importance of using a non-grey model atmosphere 
to provide the boundary conditions for evolution. 

Figure|S]shows the evolutionary tracks of models of 0.15, 
0.4, and 0.6 Mq, computed with primordial and lithium free 
element mixtures. We find that Li has little or no effect upon 
the evolution of stars of mass between 0.8 and 0.15 M©. 
The reason for this is that all the stars are hot enough to 
either ionise hydrogen or form H^, releasing more electrons 
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than can be produced by the ionisation of Li. Furthermore, 
hthium is destroyed within a few milhon years by fusion 
via the reaction ^Li(p,a)^He, which occurs before the star 
reaches the Main Sequence. As primordial stars with mass 
less than ~0.5 Mq are fully convective at zero age Main 
Sequence, all the lithium in the cooler lower mass stars is 
destroyed. 

In our earlier work, iHarris et al ] ^2004a^ we demon- 
strated the important role that has upon the evolution 
of stars of mass less than O.4M0. In order to re-asses the im- 
pact of in light of the non-grey model atmospheres we 
have computed models with and without in the equation 
of state. Figure |S] shows the evolutionary tracks of models 
of masses 0.15 and 0.4 M0 computed with and without 
in t he equation of state. Also plotted are our previous mod- 
els llHarris et al J 1200431) which are computed with a grey 
atmosphere and with the higher hydrogen mass fraction of 
X = 0.77. As the rate of the reaction H(p,e^ '^)D, is pro- 
portional to the square of the abundance of hydrogen, the 
higher mass fraction results in slightly different positions of 
the zero age main sequence on the HR diagram. Another dif- 
ference is that Harris ct al. (2004a), used the slightly higher 
value of the mixing length parameter a = 0.18. This is im- 
portant for the pre-Main and post-Main Sequence evolution, 
but has only small impact upon the Main Sequence. 

Figure |7| show evolutionary tracks computed with pri- 



mordial abundances, for 0.8 and 0.3 M© stars. One set of 
models uses ODFs, described in section 12.21 to represent 
opacity from H2, H and Li lines and the other neglects line 
opacity for these species. There is almost no difference be- 
tween the 0.8 and 0.3 M© evolutionary tracks computed 
with the H, H2 and Li ODFs and neglecting H, H2 and Li 
opacity. So, H2, H and Li lines do not have any significant 
affect on the evolutionary models presented here. 

ISiess el's! (2002) have computed zero metallicity evo- 
lution model s for the mass range 0.8 to 20 M0, and 
iMarigo et all (|200ll) for 0.7 to 100 Mq. The core hydrogen 
bur ning evolution a ry lif etimes given by ISiess et al.l (l2002h 
and iMarigQ et al.l (1200 ll) for a 0.8 Momodel are 14.8 and 
13.75 Gyr respectively. This compares with our computed 
hydrogen burning lifetime of 11.96 Gyr. The 2-3 Gyr differ- 
ence betwee n hydrogen b urnin g life times computed in t his 
work and bv lSi'ess et al] (l2002ll and iMarigo et alJ (l200l[) is 
primarily due to the use o f diff erent initial hydroge n mass 
fractions. ISiess et alJ (l2002l) and IMarigo et all i200 jl use an 
initial H mass fraction of 0.77, compared to the value of 
0.752 used in this work. 

Figure |H| shows the 14 Gyr isochrone derived from 
our models computed using primordial abundances, with 
grey and non-grey atmospheres and with and without Kg". 
The primary differences between the grey and non-grey 
isochrones occur upon the red giant branch, where mod- 
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els computed with non-grey boundary conditions result in 
lower effective temperatures for a given luminosity. Figure|U] 
shows the low mass end of the 14 Gyr isochrone computed 
with grey and non-grey model atmospheres and with and 
without Hj~ . Also shown in this figure are the isochrones of 
iHarris at alJ ll2004all and the models of SBLHB94. The addi- 
tion of H J clearly results in a significant drop in TeS the use 
of the grey model atmosphere also decreases Tcff . The mod- 
els of SBLHB94 fall between the non-grey isochrone of this 
work and the grey isochrone of Harris ct al. (2004a). Table 
121 lists broad band colours and atmospheric parameters, for 
our non-grey evolutionary models. For masses below 0.5 M©, 
the models which exclude Hg" are hotter, show bluer colours, 
and are correspondingly more luminous than those which 
account for H^. There is in general good agreement be- 
tween the colours of SBLHB94 and those computed in this 
work, although the atmospheres and colours of SBLHB94 
are slightly cooler and redder. This may be due to the fact 
that SBLHB94 use older collision induced absorption data, 
which in general gives a higher opacity and may result in 
a cooler star. The broadband color indices were compute d 
by u sing the bandpasses given bv iBessell fc Bretd Jl98d) : 
iBess ell (_199flj and were calibrated by using a spectrum of 
Vega. 



4 CONCLUSION. 

Lithium based species have been added t o our equation of 
state, enabling us to confirm the finding of lMaver fc Duschl 
i2005h that ionised lithium is the dominant positive ion at 
low temperatures, indirectly increasing opacity. This equa- 
tion of state has been included in our opacity computations, 
which indicate that the Rosseland mean opacity is a poor 
representation of opacity at low densities and for temper- 
atures below 6000 K. Our grey and non-grey model atmo- 
spheres also show that a grey atmosphere is inadequate for 
giant stars with Tcff < 6000 K. We also confirm the finding of 
SBLE[B94that the Grey approximation is invalid for dwarf 
stars with Teff<5000 K. However, we find the deviation of 
grey from non-grey model atmosphere to be smaller than re- 
ported by SBLHB94, which may be due to the use of older 
collision induced absorption data by SBLHB94. 



We have computed a set of stellar evolution models 
for primordial stars of mass 0.8-0.15 Mq, using a non-grey 
model atmosphere to provide the surface boundary condi- 
tions. The models are computed from protostellar collapse 
to the peak of the red giant branch or the point at which 
the star turns of the giant branch. These models have Tcff 
high enough (> 4000 K) to ensure sufficient electrons are 
released via the formation of H^ or H|^ , so that lithium has 
little or no effect on protostellar collapse and stellar evo- 
lution. Furthermore Li is fully destroyed before the Main 
Sequence in stars of mass less than ~0.7 M©. However, the 
role of lithium in the formation of primordial stars and low 
mass brown dwarfs, which are to cool to burn lithium, still 
needs to be investigated. 

In line with earlier work SBLHB94 we find that using a 
non-grey model atmosphere gives a lower effective tempera- 
ture and luminosity for stars of mass less than O.SMq. How- 
ever, we find the deviation of Main-Sequence grey from non- 
grey models is smaller than reported by SBLHB94, which 
again may be due to the use of older collision induced ab- 
sorption data. More importantly, our non-grey models have 
significantly lower effective temperatures for a given lumi- 
nosity, on the Hayashi track and red giant branch, than do 
our grey models. We also find that excluding Hg" from the 
computations of stellar evolution for stars with mass less 
than 0.5 M© results in cooler, bluer colours and less lumi- 
nous stars. 
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Figure 8. 14 Gyr isochrone, determined from models computed with primordial element abundances and with grey and non-grey 
atmospheres, and a grey model atmosphere with Hg" neglected from the equation of state. 
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Table 2. Broad band colours and atmospheric parameters of 14 Gyr primordial non-grey model stars, 
computed with and without Hg" in the equation of state. Also listed is the data of SBLHB94. 
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